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Abstract 

A general framework to calculate the Zener current in an indirect semiconductor with an exter- 
nally applied potential is provided. Assuming a parabolic valence and conduction band dispersion, 
the semiconductor is in equilibrium in the presence of the external field as long as the electron- 
phonon interaction is absent. The linear response to the electron-phonon interaction results in a 
non-equilibrium system. The Zener tunneling current is calculated from the number of electrons 
making the transition from valence to conduction band per unit time. A convenient expression 
based on the single particle spectral functions is provided, enabling the evaluation of the Zener tun- 
neling current under any three-dimensional potential profile. For a one dimensional potential profile 
an analytical expression is obtained for the current in a bulk semiconductor, a semiconductor under 
uniform field and a semiconductor under a non-uniform field using the WKB (WentzelKramers- 
Brillouin) approximation. The obtained results agree with the Kane result in the low field limit. 
A numerical example for abrupt p — n diodes with different doping concentrations is given, from 
which it can be seen that the uniform field model is a better approximation than the WKB model 
but a direct numerical treatment is required for low bias conditions. 



I. INTRODUCTION 

The continued device scaling in the semiconductor industry has enabled the fabrication of 
devices with nanometer sizes. The scaling of the dimensions was however not accompanied 
by a similar reduction in operating voltages and, as a consequence, large electric fields are 
present in these nanometer-sized semiconductor device structures. Accordingly, tunneling 
from valence to conduction band, also referred to as Zener tunneling, produces a significant 
current, while further reduction of the device dimensions will enlarge the role of Zener 
tunneling even more. In the end, Zener tunneling could be made so efficient that it can 
provide significant drive current for a novel type of transistor [1]. 

Zener tunneling was first described by Zener [2] and the first calculations of the Zener 
tunneling rate in a uniform electric field were made by Keldysh [3] and Kane [4] for a direct 
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semiconductor based on a two band model. In the presence of a non-uniform field however, 
the WKB approximation can be used provided the applied bias voltage is generating slowly 
varying potential profiles inside the active device area [5] . The problem of Zener tunneling in 
an indirect semiconductor was treated earlier by Keldysh [6] and Kane [7], whereas Schenk 
[8] formulated a model for the tunneling probability in a uniform field using the Kubo 
formalism[9, 10]. Another noteworthy approach is the calculation by Rivas [11] based on a 
Green's function formulation of Fermi's golden rule to calculate the current through a silicon 
tunnel diode. 

Since most of the available models are one-dimensional, the currently most popular ap- 
proach to calculate Zener current in a device is to choose a set of one-dimensional tun- 
nel paths and to determine the tunneling probability along these paths within the WKB 
approximation[12, 13]. However, as all one-dimensional models ignore the pronounced two- 
dimensional shape of a realistic potential profile, the development of a more comprehensive 
calculation of Zener tunneling in semiconductors is in order. Furthermore, Zener tunneling 
should not only be modeled accurately for two-dimensional potential profiles but should as 
well incorporate phonon-assisted tunneling phenomena. 

Various general methods for dealing with electronic transport in semiconductors taking 
quantum effects into account have been developed such as the Wigner transport equation [14- 
16], the non-equilibrium Green's function formalism [17], the Pauli master equation[18] and 
the quantum balance equations [19]. But all methods come with their own limitations and 
no comprehensive generally applicable scheme for electronic transport is available due to the 
absence of a general non-equilibrium thermodynamic framework. Furthermore all methods 
can become computationally very expensive when multiple bands have to be included as 
required by Zener tunneling. 

Being mediated by phonon scattering, the process of Zener tunneling in an indirect semi- 
conductor strongly contrasts with conventional single-band transport where phonon scat- 
tering acts as a mere dissipative mechanism limiting the drive current. This conceptual 
difference in the role of phonon scattering jeopardizes the useability of the above mentioned 
methods to model Zener tunneling[12] even further. 

In this paper we provide a framework to calculate the current due to phonon-assisted 
tunneling in indirect semiconductors. In section II, we discuss the valence and conduction 
band electrons in the indirect semiconductor under the application of an external potential. 
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In section III, we establish the Hamiltonian of the phonons and the electron-phonon inter- 
action. In section , IV we construct a non-equilibrium system starting from an equilibrium 
system consisting of three non-interacting ensembles: valence band electrons, conduction 
band electrons and free phonons, all being characterized by a grand canonical ensemble den- 
sity matrix of the Gibbs form. We use perturbation theory to calculate the non-equilibrium 
density matrix up to first order in the electron - phonon interaction, the latter mediating 
the band-to-band tunneling events. Then, we extract an expression for the steady state 
phonon-assisted tunneling current from the non-equilibrium density matrix as well as the 
net rate of electrons making the transition from valence to conduction band. In particular, 
for the sake of computational efficiency, we rewrite the tunneling current in terms of properly 
chosen spectral density functions. Using the proposed method, an analytical expression is 
derived for tunneling in the presence of a one-dimensional potential profile in three cases: 
1) no external electric field, 2) a uniform electric field and 3) a non-uniform ID field within 
the WKB approximation. Finally, we discuss both the analytical results and the numerical 
calculations. 

In a recent paper [20] we apply the presented formalism to the case of a tunnel field-effect 
transistor and show that the use of semi-classical model introduces very large errors due to 
quantum confinement near the interface. 



II. ENVELOPE FUNCTION APPROXIMATION FOR AN INDIRECT SEMI- 
CONDUCTOR 

In this section we consider the valence and conduction band electrons in an indirect semi- 
conductor subjected to an external bias voltage. We introduce the single-electron energies 
and eigenfunctions required to construct the non-interacting electron Hamiltonian and some 
related second-quantized operators. In section IV, the obtained operators are used to con- 
struct the non-equilibrium density matrix by introducing the interaction with the phonons. 

In a semiconductor an electron is subjected to a periodic potential energy £/i at (r) caused 
by the ions and electrons in the lattice. In a semiconductor device, doping and bias voltages 
externally applied at the contacts induce an additional non-periodic potential energy f/ ext (r). 
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The single-electron Schrodinger equation reads 

h 2 v 2 



2m 



+ C/ lat (r) + C/ ext (r)U(r) = ^(r). (1) 



In the absence of the non-periodic potential, the periodicity of the potential results in 
Bloch wavefunctions 

h 2 v 2 \ 

+ ^iat(r) exp(ik • r)u nk (r) = £„ k exp(ik • r)u nk (r). (2) 



2m 

where n is a band index, k is a wavevector and E nk describes the energy dispersion relation. 
For semiconductor device applications, n is essentially running over c and v, respectively 
denoting the conduction and valence band separated by the bandgap, while the states that 
participate in the conduction process are located near the Fermi level. In a degenerate 
semiconductor, as often employed when studying Zener tunneling, the Fermi level will be 
located inside the conduction or valence band. 

In the indirect semiconductor under consideration, we assume a single valence band max- 
imum occurring at k = and one of the conduction band minima being located at k = k , 
as shown in Fig. 1. The energy dispersion can be approximated around the conduction band 
minimum and the valence band maximum using the effective mass approximation: 

h 2 

E vk « E v0 - — r|k| 2 fork^O, (3) 

2m* 

h 2 

E ck « E cko + — — |k - k | 2 for k « k . (4) 

2m* 

When an external potential is added, the wavefunctions are no longer Bloch functions 
and no energy dispersion can be defined. However, the Bloch states solving Eq. (2) still 
provide a complete orthogonal set of basis states in which each one-electron eigenstate ^(r) 
can be expanded as 

V>(r) = ^ a «k exp(ik • r)u nk (r) (5) 

n,k 

On the other hand, we may as well expand ^(r) around a high-symmetry point within the 
first Brillouin zone, such as the T-point (k = 0), writing 

^( r ) = ^2xno(r)u n o(r) (6) 

n 

where Eq. (6) defines the envelope functions Xno( r )- 
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FIG. 1. The Si bandstructure along the [100] direction[21] around the valence band maximum and 
conduction band minimum indicating the minimum of the conduction band at k = ko 

For the conduction band electrons, the envelope functions denned in Eq. (6) are rapidly 
oscillating since the conduction band minimum has a large wavevector. In order to end up 
with smoothly varying envelope functions, we may alternatively define the set Xnk by 



^(r) = Yl Xnk„(r) exp(ik • r)w nko (r). 



(7) 



For a given eigen energy E, the dominant contribution in Eq. (7) comes from the component 
with its eigen energy E nko closest to E. By restricting the summation of the index to the 
conduction band, the wavefunction of a conduction band electron can be written as 



^c(r) = Xc(r) exp(ik • r)w cko (r) 



(8) 



with a slowly varying envelope function, i. e. the fourier components a„k are only significant 
for k k . Since we consider the conduction band to be a non-degenerate valley at k , we 
conclude that -u c k ~ w c k for k ps k and 



Xc(r) « ^ fl ck exp (i(k - k )-r) 



(9) 



Moreover, adopting the effective mass approximation, we may approximate the kinetic plus 
the periodic term of Eq. (1) using Eq. (2) to get: 

h 2 v 2 

- + fiatlrj j Xc{r) exp [ik • r) u cko {r) = ^ 

k 



2m 



+ Ui at (r) J Xc(r) exp (ik • r) u c]l0 (r) = ^ a c k^ckexp (i(k - k )-r) exp (ik • r) M c k (r) 



h 2 v 2 

2m* 



Xc(r) exp (ik • r) w cko (r) (10) 
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The envelope functions are determined by combining Eq. (1), Eq. (8) and Eq. (10), which 
results in the effective Schrodinger equation 

(Eco ~ ^V 2 + C/ ext (r)) Xc (r) =E Xs (r). (11) 

For valence band states the interaction between the different degenerate valence bands 
must be taken into account. Nevertheless, we will restrict ourselves to the effective mass 
approximation also for the valence band which is a common approach [8]. To go to a more 
realistic valence band structure, 6-band kp theory could be employed without compromising 
the following approach to calculate the current. 

A major improvement of going beyond the effective mass approximation will be a more 
accurate description of the complex band structure in the bandgap which affects the decay of 
the wavefunctions penetrating the bandgap. Inspection of the complex band structure for the 
case of silicon reveals that the effective mass approximation overestimates the decay. Hence, 
the real tunneling probabilities will be higher than those predicted by the effective mass 
approximation. A second consequence of taking non-parabolicity into account is a different 
density of states, which will greatly influence the valence band because of its anisotropy and 
the presence of the split-off band. 

In the effective mass approximation the wave functions for valence band electrons are 

?/v(r) = Xv(i> v0 (r) (12) 

where the envelope functions are determined from 

(e v0 + ^V 2 + C/ e xt(r))xv(r) = E Xw {v). (13) 

We can now construct the electron basis set consisting of the valence and conduction 
band eigenstates 

Mr) = Xve{r)u v0 {r) (14) 
0rf(r) = Xc£(rK k() (r) exp(ik • r), (15) 

where XvK r ) and Xd{ r ) are f ne envelope eigenfunctions of Eq. (11) and Eq. (13) with 
eigenvalues E vi and E C £ respectively. £ is a subband index running through the set of 
relevant quantum numbers. 
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The one-electron states allow us to write the electron Hamiltonian in second quantization 

H c \ = ^2 E vetyb e + E c£ c\c e (16) 
t 

where bg and q annihilate an electron in the states v ^(r) and (f)a( r ) respectively. Further- 
more, it proves convenient to define electron field operators 

Vv(r) = ^<Mr)6,, (17) 

e 

&(r) = J2M*)ce. (18) 

i 

to represent the electron density and other position dependent quantities. As explained in 
the next section, the tunneling current is related to the expectation value of the valence and 
conduction band electron number operators 

N v = ^b\b e = I d 3 r^(r)^(r), (19) 
N c = J24c £ = I d 3 r^(r)4,(r). (20) 

III. PHONONS 

In general, one may want to consider various agents other than phonons that are proven 
to assist Zener tunneling, such as traps and defects, but in this work we have investigated 
the role of phonons only. 

The free phonons are described by the Hamiltonian 

#PH = ^ ^qa^Oq (21) 

q 

where a q annihilates a phonon of mode q with energy hw^. 

As this work focuses on indirect semiconductors, we specifically consider the interaction 
between electrons and short-wavelength phonons bridging the gap between the top of the 
valence band and the bottom of the conduction band located at k . The corresponding 
interaction hamiltonian reads [22] 

**' = Yl (?vtcej%c e (a q + ot q ) + h.c.) . (22) 
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The coupling strengths g v fc£'q incorporate all possible interband transitions and typically 
reflect the features of the envelope wave functions. Alternative interaction mechanisms ca- 
pable of assisting interband transitions such as localized trapped electrons are not considered 
in this work. 

IV. NON-EQUILIBRIUM DENSITY MATRIX 

In this section we construct the non-equilibrium density matrix by tracing the time evo- 
lution of the density matrix up to the lowest order in the electron-phonon interaction. Me- 
diating electron transitions between the valence and conduction band, the non-equilibrium 
density matrix carries a non-zero tunneling current. 

A. Equilibrium 

In the absence of interactions we may consider the valence band and conduction band 
electrons as well as the phonons three separate entities that are individually characterized 
by their equilibrium density matrices. The density matrix governing the entire, uncoupled 
system therefore reduces to the direct product of the three density matrices, which can be 
written in its grand-canonical form 

Qo = Qv <8> Qc <8> £?ph (23) 

exp (-p(H A - /i v iV v - fi c N c + H PH )) 
= * % L (24) 

where /3 = l/(kT) and Z denotes the partition function to ensure that Tr(£ ) = 1- The 
chemical potentials of the conduction and valence band electrons fi v ,fi c will be determined 
by imposing charge neutrality in the device contacts whereas the phonon system has zero 
chemical potential. 

Accounting for spin the net charge density of the equilibrium system can be calculated 
from 

p net (r) = -2eTr ((^(r)4(r) - ^ v (r)^(r))^ ) ■ (25) 
In order to facilitate both numerical or analytical evaluations, we define the valence and 
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conduction band spectral functions based on the envelope functions: 

A v , c (r,r';E) = 2ir6(E-H Y , c ) (26) 
= 2tt XvA*)5(E ~ E v ,a)x:, c e(r') • (27) 

£ 

The net charge density averaged over a unit cell can now be evaluated as 

p net (r) =-2e J ^(f c (E)A c (r,r-,E)+(l-f Y (E))A v (r,r-,E)y (28) 

where the valence and the conduction band distribution function is a Fermi-Dirac distribu- 
tion 

fvAE) = exp(/3(£-^ c )) + l- (29) 
B. The quantum Liouville equation 

Taking H to be the Hamiltonian of the non-interacting electrons and phonons, i.e. 

H = H A + H PU , (30) 
we express the time dependent density matrix in the interaction picture as 

Q(t)= exp(iH t/h)g(t)exp(-iH t/h) (31) 
to calculate its time evolution from 

ih^g(t) = [H'(t),g(t)} (32) 
where H'(t) denotes the electron - phonon interaction, expressed in the interaction picture, 

H'{t)= exp(iH t/h)H'exp(-iHot/h). (33) 
Integrating Eq. (32) and keeping only terms of first order in H'{t), we obtain 



l 



Q (1 \t) = g -- dr[H\r),g ] (34) 

since £> commutes with H and is unchanged in the interaction picture. 

Since H' represents a phonon creation or annihilation event and £o is diagonal in the 
phonon states, the diagonal elements of the density matrix are not affected to first order. 
Hence, the electron density remains unaltered and can still be calculated from Eq. (28). 
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V. THE STEADY STATE CURRENT 



The quantity we want to calculate is the current I carried by the interacting electron - 
phonon system at steady state, i.e. when t — > oo. As a definition, the tunneling current 
basically represents the rate of change in the number of electrons in the conduction band 
per unit time: 

I v = -2e lim — Tr (N v g(t)) , (35) 
t^oo at V / 

I c = -2e lim A Tr (N c g(t)) , (36) 
t^oo at \ / 

where I v + I c = holds due to particle conservation and the factor of 2 is due to spin. 

Since the valence band number operator N v commutes with H , N v does not gain any time 
dependence in the interaction picture. Making use of the Liouville equation and Eq. (35), 
we obtain 

J v = -^Hm jf drTr^N v [H'(t),[H'(T),g }]y (37) 
Exploiting the identity 

Tr(A[B, [C, £>]]) = Tr([A, [B, C]]D) (38) 
which follows from the invariance under cyclic permutations, we rewrite Eq. (37), 

I y = | hm ^drTr {[[N v , H'(t)\, H'(r)] Qo ) . (39) 
A. The phonon-assisted tunneling current 

For the phonons, the interaction Hamiltonian in the interaction picture is given by 
H'(t) = ^2 9vkd>'$th> (aq exp(i(£ v£ - E ci , - hu^t/h) 

+ al q exp(i(£ v/ - E cfJ + fkjjjt/hj) + h.c. (40) 



under the assumption cu_ q = w q . 

The commutator [N v , H'(t)] can be evaluated straightforwardly to yield 

[N y ,H'(t)] = ^ g v ece>qfyce> (a q exp(i(E v£ - E ci > - hu^t/h) 



+ a f _ q exp(i(£ v£ - E c& + fiw^t/h)) - h.c. (41) 
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Since go is diagonal in the phonon operators, only terms connecting a phonon annihilation 
operator in [N v , H'(t)] with a creation operator in H'(r) and vice versa do contribute to the 
average of the double commutator appearing in Eq. (39). The resulting expression for the 
current therefore simplifies to 



7 V =-| lim / dr Tr ( 

J U 00' r. 



2COS ((E V £ - E cV - flWq)(t - t)/K) Qq 



+ 



SWq&{Q'aL q , (gvtct'JiceaL^ 2cos ((E ve - E cSJ + hu)q)(t - t)H) g^j . (42) 



Bearing in mind that the initial equilibrium state describes uncoupled electrons and phonons, 
we may employ Eq. (23) to simplify the evaluation of the traces. For example, 

Tr (b\c e a q a q ,cJ'&rf?o) = Tr (b\b r g^j Tr (q/cJ,£ c ) Tr (a q a q/ £ PH ) 

= M*'r'<W/v(#v<) (l ~ /c(^')) ( X + K^q)) , ( 43 ) 

where v(E) denotes the Bose-Einstein distribution 

1 



u(E) = - 
K ' exp(/3£)-l 

Taking the limit t — > oo and using the identity 



(44) 



lim / dr cos cut = ir5(u), 



(45) 



we derive the steady-state current 
47re 



J v — r~ ^2 \9v£c£'q\ 2 

X (/ v (Erf)(l - / c (£rfO)M^q) + 1)^(^ " " ^q) 

- f c (E ce r)(l - f v (E ye ))l/(fiWq)5(E ve - Ed' - ftlUq) 

+ f v (E ve ){l - UE^vihu^iEvt - E ce , + fiwj 

- f c (E cf ,)(l - / v (^))(i/(7ia; q ) + l)<5(S v/ - E cl , + 7ia; q )). (46) 

The four terms contributing to the current in Eq. (46) can be interpreted in a similar way as 
in Ref. [11]. For instance, the first term refers to an electron being excited from the valence 
to the conduction band while emitting a phonon, etc. 
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B. Evaluation of the current using spectral functions 



Alternatively, the tunneling current may be written as an energy integral 

i = -| / ^ E ((/v(£)(i - ME - ^ q ))(^ q ) + 1) 

q 

- / c (£ - «w q )(l - / v (£)M^ q ))T v cm (£,q) 

+ (/ v (£)(l-/ c (£ + ^ q )M^ q ) 

- / c (£ + /^ q )(l - / ¥ (S))(i/(fc; q ) + l))T v abs (£, q)) (47) 

where the property of the delta function 

/+oo 
dE 5{E - E cl > - hu^) 5{E - E v£ ) 
■oo 

has been used and the probability of exciting an electron from the valence band to the 
conduction band with emission or absorption of a phonon has been introduced: 

T abs,em ( ^ q) = ^2 ^ \ g ^,^^ E - E v£ ) X 5(E - E cf , ± fiwj. (48) 

w 

Equally it would be possible to write Eq. (47) in terms of the reverse transition by using 

Tr(E, q) = T v abs (£ - hu v -q), (49) 
T abs (£, q) = Tr(E + /^ q , -q). (50) 

Further numerical processing of the current formula requires an explicit representation 
of the scattering matrix elements. The latter are related to the bulk coupling strength M q 
through 

#vfc£' q = M q J d 3 r (j>l e (r) exp(iq • r) c ^(r). (51) 
In this work, we have borrowed M q from the deformation potential interaction [22] 



M ^ D ^2^i < 52) 

where p s , D and Vt respectively represent the the semiconductor mass density, the deforma- 
tion potential and the total volume. In order to express the matrix elements in terms of the 
envelope functions, we assume that the basis functions are normalized on the lattice unit 
cell fl c and define 

M q = M q [ d 3 r< (r) Mck() (r). (53) 
13 



Since the envelope functions are slowly varying in space, they can be taken constant along 
each unit cell and the matrix element from Eq. (51) can be approximated as 

Wq = Kj d3r ^( r ) exp(i(q + k )-r) X d'(r). (54) 

Again, within the perspective of numerical evaluation, we may use the single-particle 
spectral functions A(r, r'; E) defined in Eq. (27) to rewrite the tunneling probability 

T v abs ' cm (£,q) = |M;| 2 Jd 3 rjd 3 r'(exp(i(q + k )ir' - r)) 

xA v (r,r';E)A c (r',r;E±hujS)- (55) 

As the matrix elements g v fcfq turn out to be strongly peaked for q = — ko, we can assume 
negligible dispersion for the phonons assisting the tunneling processes, i.e. w q ~ w_k = O0k - 
Hence, the occupation probabilities in Eq. (47) become independent of q and we can define 
effective transition probabilities by 

T v abs ' cm (£) = T v abs ' cm (£, q), (56) 
q 

to write the tunneling current as follows: 

/ = ~f / ^((/v(^)(l-/c(^-^ ko ))K^ k o) + l) 

- f c (E - hw ka )(l - f v (E))v(hw k0 ))Tr(E) 
+ (f v (E)(l-f c (E + tu^))v(tiw k0 ) 

- f c (E + fc^l - f v (E))(u(Huj k0 ) + l))T v abs (£)Y (57) 



Furthermore, neglecting the dispersion of in Eq. (55) for analogous reasons, we may 
carry out the sum over the wave vectors, yielding 

^ exp(i(k - q) • (r' - r)) -+ nj ^exp(i(k - q) • (r' - r)) = fi5(r - r'). (58) 
In turn, the tunneling probability of Eq. (56) simplifies to 

T abs,cm (jE) = n \ M ^fJ d 3 r Av(r? r; E) A ^ r . E ± ^J. (59) 

involving only the diagonal terms of the spectral density. 
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VI. A ONE-DIMENSIONAL POTENTIAL PROFILE 



To demonstrate the outlined method, we consider the case where the external voltage 
gives rise to a one-dimensional potential profile, say in the x-direction, while translational 
invariance exists in the y- and z-directions. 

The valence and conduction band envelope Schrodinger equations read 

h 2 d 2 



2m v dx 2 



Ejr(K) + U ext (x)) X v(r;K,E) 



= E Xv (x;K,E), (60) 
(E g -^ + E-(K) + U^))x c (r-,K,E) 

= E Xc (x;K,E). (61) 

Here, the valence band maximum and the conduction band minimum respectively correspond 
to the energy values E = and E = E g , whereas K denotes the two-dimensional wave vector 
reciprocal to R = (y,z). Representing the kinetic energies in the transverse directions, 
E^ C (K) can be expressed in the effective mass approximation as 

^(K) = (62) 

Based on the solutions of the Schrodinger equation, the diagonal part of the spectral densities 
read 

4r,c(r, r; E) = / --| Xv , c (r; K, £)| 2 (63) 
provided that the wavefunctions are delta normalized 

/+oo 
d 3 r X ; c (r; K, £) Xv , cc (r; K', E>) = (2tt) 3 5(E - E')6(K - K'). (64) 
-oo 

Changing integration variables in Eq. (63) from K to E^ c and the polar angle of K, we 
obtain 

m f°° rl F 1- *- 

A ViC (r,r;E) = ^y o ^HXv,c(r; E^, E)\ 2 (65) 

assuming |Xv,c( r ; K, E)\ 2 only depends on the magnitude of K. Note that we have assumed 
an isotropic mass. In case the mass is anisotropic with a mass m v>x , ra V;2/ , m V;2 in the x, y 
and z direction respectively, m v in Eq. (65) has to be replaced with ^m v ^m v ^ z while m v in 
Eq. (60) has to be replaced with m V:X . 
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A. Bulk 



As a first example, consider a bulk semiconductor where no external field is applied, 

£4xt(r) = 0. (66) 

Clearly, no Zener tunneling is expected in the absence of an electric field but the treatment 
of an unbiased bulk semiconductor illustrates the consistency of the formalism. 

The normalized solutions to the envelope Schrodinger equation for E < E^ (K) are 

±ixJ '-2m v /h 2 (E+E^-(K)) 

X±(r; K, E) = 6 ^e iK - R , (67) 

V-ff^ + ^K)) 

±ix^/2m c /^(E-E g -E^CK)) 

Xt(r; K, E) = - ^e iKR (68) 

^(E-E g -Ei-(K)) 

where the superscript ± denoting right and left running waves had to be added to complete 
the basis. The valence band spectral function for E < is 

A v (r, r; E) = ^jf ^(| x +(r; E£, E)\ 2 + | x ;(r; E v \ E)f) 



m v / m 



E 

dEd 



h 2 2ir V 2h 2 Jo v \/—{E + E^) 

The conduction band spectral function for E > E g is obtained by substituting v with c and 
E with (E g — E) in the envelope functions 



trie 



A c (r,r;E) = ^- x /2(E-E g ). (70) 



The transition probability for E g — fku < E < is 



TT(E) = 0, (71) 

3 

,2(m v m c )2 



T* hs (E) = tf\M^ r^y- y/-E(E + nw- E g ). (72) 

Eq. (71) reflects that no transitions from valence band to the conduction band by phonon 
emission are possible in bulk semiconductors (E g > 0). 

In Eq. (72) the condition for transitions to occur is that the phonon energy exceeds the 
bandgap which is generally not the case in semiconductors. As expected, no phonon-assisted 
tunneling is present in the absence of an external potential. 
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B. The ID Uniform field 



A situation more relevant to the problem of Zener tunneling is the application of a uniform 
electric field or a uniform force F to an indirect semiconductor 

U cxt (x) = -Fx. (73) 

The solutions to the effective mass Schrodinger equations are written in terms of the Airy 
function defined as the bounded solution of the differential equation [23] 

d 2 

— — Ai(x) = xAi(x) (74) 
ax 2 

and Ai(0) = 3- 2 / 3 /r(2/3). 

The delta normalized solutions to the valence band Schrodinger equation are 



ME)= ^ Ai (± + ^m ]e ^ (78) 



2vr / x E — Eg — E^~ (K) \ iK-R 



Xc(r;K,E) = ,I—Ml-- ) c (76) 



with 

^ = s£f < 77 » 

and where the orthogonality of the Airy functions is taken from Aspnes[24] 

du Ai(x + u)Ai(y + u) = S(x - y). (78) 

' — oo 

The spectral functions are 

my I f 00 , > .2 fx , £ 



J — c 



- 4 -(" :£ »=#^X d " Ai U + ^ + "J- < 79 > 

^■^'=^r d " Ai2 (-i + ^? + ")' <8o) 

Since the field is uniform, the transition probability is independent of the energy 
T* hs ' em (E) = T^ hs ' em (0). The integration of the spectral functions yielding the transition 
probability is straightforward using Eq. (B34b) from Aspnes[24] and given in the appendix, 
the result is 

ir-m = w A ^: F j A k ( ^mM } (81) 
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with A the device area in y and z direction, 1/m— l/m v + l/m c and Ai 3 (x) is defined as 

Ai 3 (ar) = ]- (Ai(x) + xAi'(x) + a; 2 Aii(a;)) (82) 

with 

/■oo 

Aii (a;) = du Ai(u). (83) 



The transition probability can be approximated for large fields using Eq. (9) from Heigl 
et al.[25] 

Ai 3 (x) w — ^exp — 12 (84) 



2^ 3 
resulting in the Kane expression: 

r - nB) „ W^zlM exp 4 VM^Ml (85) 

C. ID Non-uniform fields: the WKB approximation 

Next, we calculate the phonon-assisted tunneling current in a one-dimensional device 
structure with external potential 

U cxt (r) = U cxt (x) (86) 

using the WKB approximation [26]. The envelope functions within the WKB approximation 
are 



exp (- ( x dx'J^(E - U ext (x>) + E±)) 
Mr; K, E) = ^ J - 2 V ^ Y V-», (87) 

where x t v,c indicate the location of the turning points as illustrated in Fig. 2. In the spirit 
of the WKB approximation, the normalization can be verified by comparing to the uniform 
field case. The turning points are determined by the equations 

E - U cxt {x tv ) + Ey = 0, 
E g + U cxt (x tc ) — E + E^r = 0. (89) 
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FIG. 2. Illustration of the turning points x tv and x tc for E = appearing as the solutions to 
U ext (x tv ) — = and E g + U cxt {x tc ) + E^ = respectively. 

In appendix B, an approximate expression for the spectral functions and the tunneling 
probability is obtained. The results are 

A (r r- E) ~ <=x P (-2 d,- K .(,'; £)) 



exp(2r :r dx'/tcfa;': i?) ) 



and 



Afi|M4 o | 2 (m v m c )2 



2 4 7T2^2(£' zp hw) 41714 y/Ul 



ext V"^max 



± 



X 



ex p (- 2 L:r d ^c(^;^)) 
r:r da;/ /^(a:';i?) 



ex p ( 2 ix t zr dx ' k A x '; e ± ) 

f^™ ax dx'/ k c (x';E± hw k() ) 92 
where x tV)C o are the turning points x tVjC for E^ c = 0. x = x^ ax denotes the point where the 
imaginary wave vector of the valence and conduction band equal each other: ^(x^; E) = 
K c(^max5 E ^ ^ko)- ^cxt (^max) denotes the first derivative of the external potential at x — 
x max- The imaginary wavevectors defined as 



'2m 

Kv (x;E) = \l-^(E-U„ t (x)), (93) 



1 2m 

k c (x; E) = \J-f^{E g + U ext (x) - E). (94) 

Taking the limit of the tunneling probability for small uniform fields U cxt (x) = —Fx, the 
WKB expression approaches the uniform field expression obtained in the previous section. 
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D. Non-uniform fields: numerical evaluation 



In a real semiconductor under a one dimensional external potential, the external potential 
can be taken uniform in two contact regions on the right and left hand side. The discretized 
wavefunction can be determined by discretizing the Schrodinger equation and applying so- 
called transmitting boundary conditions [27]. In a Zener tunnel junction, the wave will not 
be transmitted but completely reflected by the potential barrier. The wavefunction will only 
have a decaying component in the contact opposite to the side where the wave is injected. 

The incoming component of the wavefunction is given by the bulk wavefunction from 
Eq. (67). The spectral functions are obtained by integrating the wavefunctions with respect 
to energy according to Eq. (63). 



VII. DISCUSSION 
A. Generation rate 

Rather than the transition probability T^ hs,em used in this paper, the so-called Band- 
to-Band Tunneling (BTBT) generation/recombination rate per unit volume G is usually 
presented to calculate the Zener tunneling rate. The current as a function of generation rate 
is given by 

I = I™ + J v abs (95) 

with 

jabs.cm = - e J d 3 r Gf s ' em (r) 

(/v(C/ ex t(r))(l - / c (C/ext(r) ± hwv«)){\ t\ + K^kJ) 

- / c (C/ext(r) ± ^ ko )(l - / v (C/ ext (r)))((^ ± i + 1/(710*,)). (96) 

G abs,em (r) denotes the generation rate for an electron with a tunnel path starting at r and 
ending at r + l abs > cm ( r ) such that C/ ex t(r) = E g =f hu ko + C/ ext (r + ^^(r)). 

Comparing the semi-classical equation with the expression for the current, the relation 
between the transition probability and the generation rate is given by 

G v (x) = l l^M T ^(U ext (x)). (97) 
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In the low field limit (Kane/Keldysh result) and assuming an anisotropic effective mass, 

(98) 



^abs,em ^abs,em p 5 [ 4 ygmgE^Mf 
G v,Kanc = G b 2 eX P I ~g ^ 



with 



G abs,c m = ^ | M / o 1 2 j ^W^^ (9Q) 



where (70 is a degeneracy factor. The factor of 2 for spin degeneracy has to be accounted 
for in g together with the appropriate valley degeneracy and the phonon degeneracy. The 
generation rate in the uniform field case is 

G abs, Cm = G abs,c mF | /2 hmi(E^Thw) ^ (l0Q) 



(hF) 



3 



with 

h(x) = 2^x*Ai 3 (x). (101) 

In Fig. 3, the generation rate is plotted as a function of the electric field for Si using the 
uniform field model and the low field limit. The degeneracy factor is 2 due to spin, 4 due 
to the four X valleys with their transversal mass in the [100] direction and 2 due to TO 
phonon degeneracy. Since the tunneling probability is strongly dependent on the effective 
mass the tunneling to the two X valleys with their longitudinal mass in the [100] direction 
can be neglected. Similarly the tunneling from the heavy hole band can be neglected. At 
high fields, the error of the Kane model due to the use of the asymptotic expression for 
Ais(x) can be observed. For the maximum field used in the plot (lOMV/cm), the difference 
between both models is a factor of 1.86. 



B. Comparing the different models for a p — n diode 

In an abrupt p — n diode, the Shockley approximation can be used to determine the 
potential profile consisting of two parabolic sections as shown in Fig. 4. The fermi levels are 
determined by imposing charge neutrality in the contacts where the carrier concentration 
can be calculated from the bulk spectral functions given in Eq. (79) and Eq. (80). 

A comparison between the numerical results and the uniform field model, the limit of 
the uniform field model and the non-uniform model using the WKB approximation is given 
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Electric field (V/cm) 



FIG. 3. Generation rate calculated using the uniform field model (solid line) and the Kane limit 
(dashed line) for Si with E g = 1.12 eV, m v = (.16, .49, .49)m and m c = (.1905, .9163, .9163)m for 
the electronic parameters and .Dko = 6 x 10 8 eV/cm, Huj = 57.6 x 1CP 3 eV, psi = 2.328 cm -3 for 
the phonon interaction [11]. Assuming tunneling along [100], the total degeneracy prefactor is 16. 
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f U ext (x) - 









' 5 2 4 6 8 10 12 



x (nm) 

FIG. 4. Parabolic potential as obtained in a abrupt p — n diode using the Shockley approximation. 
Donor and acceptor doping levels are 10 20 cm -3 , the dashed line indicates the junction. 



in Fig. 5. For the uniform field model in the presence of a non- uniform field, the usual 
approximation 



F w E g /l tnn (102) 
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with U cxt (x) = U cxt (x + hun(x)) + E g was used. The non-uniform model from this paper is 

^abs,em ^»abs,em ^ (-^g ~F f^\e.o} m x 

U v,NU - ^0 Iff 7 r 

m v,x m c,x V t^ext l^maxj 



X 



exp (-2 /* ™ dx'K c (x'; £)) 



exp (2 r^™ ax dx'/t v (x'; £ ± /^fc ) / 

r;-dxVK C (x';E±^ k0 ) (103) 

The Schenk model [8, 25] has a different dependence on the field strength (F 9//2 ) compared 
to the Kane/Keldysh/WKB result and for this reason we will not compare with the Schenk 
model. 

In Fig. 5 it can be seen that of the three semi-classical models the uniform field model 
(G Vi u) is the best. Despite disregarding the non-uniformities, the uniform field model is 
a better approximation than the WKB model and this can be explained as follows. At 
low bias (Vp n ~ 0), taking the integration limit |K| — > oo is too optimistic resulting in an 
overestimation of the current in all semi-classical models. But on the other hand the uniform 
field model is pessimistic about the tunnel probability compared to the WKB approximation 
and this partly compensates for the overestimation of the integration limits. In the high bias 
regime (V pn ^> 0), the approximation of the potential relevant for tunneling by a uniform 
field is reasonable since the field does not change much in the region where the field is the 
highest. This can be verified by the similarity of the Kane model and the WKB models. 
The difference between the uniform field model and the Kane model is that the integrals 
over x and K are treated exactly in the uniform field case based on the properties of the 
Airy functions. 

Nevertheless, for the degenerately doped semiconductors, none of the semiclassical models 
are a good approximation in the forward bias regime. When an accurate estimation of the 
current magnitude or shape for the forward bias or the low bias regime is required, a direct 
numerical procedure is required. 



VIII. CONCLUSION 



A general formalism to calculate the phonon-assisted Zener tunneling current in non- 
uniform fields using spectral functions was presented. The expression using the spectral 
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FIG. 5. Calculation of the current for a symmetrically doped Si diode with a doping level of 
5 x 10 19 cm~ 3 , 10 20 cm -3 and 1.5 x 10 20 cm~ 3 . The same Si parameters were used as in Fig. 3. 
The different current curves are obtained using the semi-classical Kane (Eq. 98, solid line with +), 
Uniform (Eq. (100), solid line with o) and WKB (Eq. (103), dashed line) model with respect to 
the direct numerical calculation of Eq. (59) (solid line). In the case of the highest doping level the 
Kane and WKB model coincide. 

function enables the calculation of the tunneling probability to two- or three-dimensional 
potential profiles. The evaluation of the current can be done using the expression for the 
current (Eq. (57)) and the expression for the transition probability (Eq. (59)). 

In the presence of a one- dimensional external potential, an analytical expression for the 
tunneling can be obtained for uniform fields using Airy functions and non-uniform fields 
using the WKB approximation. It is shown that in the weak and uniform field limit, the 
formalism from this paper reduces to the indirect Kane result. Furthermore, an improved 
uniform field model is derived and an expression using the WKB approximation are obtained. 
Comparing the different approximations in the case of a p — n junction, the improved uni- 
form field model is shown to be a better approximation than the model using the WKB 
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approximation. But for low bias conditions, no approximate model accurately describes the 
tunneling current and a direct numerical treatment is essential even in the one-dimensional 

case. 
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APPENDIX A: TUNNELING PROBABILITY IN A UNIFORM FIELD 

In this appendix, we derive the transition probability for the uniform field and a non- 
uniform field starting from Eq. (59). Using the translation invariance in y and z direction, 
the integration in y and z corresponds to a simple multiplication by the area A: 

/oo 
d^ v (r,r;0)A c (r,r;±M 
-oo 

Introducing the spectral functions for the uniform field from Eq. (79) and Eq. (80) the 
tunneling probability reads 



T^(E) =n|<| J ^ J™ dx J™ d Ul J™ duiAfix/xv+uJAi 2 (-{x+ ±hUJ ^ F S /^ + u^ 

(.104) 

The integral over x can be performed using Eq. (B34b) from Aspnes[24] 

„m m A f°° f°° x ( 2hx v /x c u 1 - ±r,UJ ^o~ E s + \ 

r v abs -™(£) = n|M;j 2 ^ duj du 2 — %=m 1 i 1 ' Fxc 



h 4 x v x c J J 4tt^/x v /x c \ (1 + {x v /x c f) 

, , m c m v A [°° f°° A f2l(x vUl + x c u 2 + ^p^[ 



3 



^v^c Jo Jo 



= n\ML t \ \ ; / dui / d M2 Aix 



(.105) 
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Defining the integral 

/»oo roc 

Ai 3 (x) = / d-ui / du 2 Ai 1 (x + ui + u 2 ) 
Jo Jo 

[~ d(u 1 + u 2 ) d(m- M2 ) A . 

= / / 7k Ali(x + Ui+U2) 

JO J-(«i+« 2 ) 

dMwAii(x + u) 



poo 

Jo 





= - (Ai(x) + xAi'(x) + x 2 Aii(a;)) (.106) 
where the last equality comes from Schenk[8], the transition probability can be written as 

3 , „.3> 2 



rpabs,em/ ztt\ Olfl/f' 1 2 ra c 771 v A (aJ v + aJ c )3 



Substituting the values for a; c and x v , Eq. (.107) reads 

3 ,-, 1 



^i ^iy i. MOT) 



7-—(B) = niKI' ^^I^ Ai, ( | . (.108) 



8-Kh 4 h3m3 \ Fsh 



3 



with 1/m — l/m v + l/m c . 

APPENDIX B: TUNNELING PROBABILITY IN A NON-UNIFORM FIELD US- 
ING THE WKB APPROXIMATION 

In the WKB approximation, the integral over E^ c yielding the spectral functions is gen- 
erally a nonelementary integral. To obtain a simplified expression, approximate the spectral 
functions by making the first order Taylor expansion of the argument of the exponential in 
Eq. (87) at = 0, 



j* dx' ^^(E-U ext (x') + E^) 
~ J* dx'^^(E-U cxt ^')) 



e' r . [2^ i 

+ -f / dx\ _. v ■= .109 
2 J X1 V ft y/E-U^ix 1 ) 



Introduce the imaginary wave vector 



2m 

k v (x;E) = X l-^(E-U cxt (x)), (.110) 
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and take the variation of the exponential to be much faster than that of the denominator in 
Eq. (87). The spectral function is approximated as 

^ r ' r '^ * 4nh^ v (x;E)f: tvo dx>/K v (x>;E) 

where x tVjC0 are the turning points x tVjC for E^ c = 0. The conduction band spectral function 
is given by 

A (r r E) ~ ^ exp(2 J^dx'^x'; E)) 

** r ' r '^ ~ Anh^ c (x;E)j: tco dx'/K c (x';E) ^ 

with 



k c {x- E) = yj -^(E g + U cxt (x) - E). (.113) 

The product of the spectral functions appearing under the integral in the expression for 
the tunneling probability 

/oo 
dxA v (r,r;0)A c (r,r;±nw) (.114) 
-oo 

can be rewritten as a prefactor f(x) and an exponential e 9 ^ 

A v (r,r;0)A c (r,r;±M = f(x)e 9 ^. (.115) 
Observing the exponential terms in Eq. (.111) and Eq. (.112), 

l>X PX 

g(x) = 2 / dx'K c (x'; E ± hou ko ) -2 dx' K v {x E). (.116) 

J Xtc J Xtv 

The exponential will decay rapidly from its maximum which is found at x — x max , which 
is determined by the condition = 0: 



9fh 



The argument of the exponential can be expanded to second order at x T 

P 3? max f^max 

g(x) « 2 / dx' k c (x' ; E ± huk ) — 2 / dx'K v (x';E) 



X 2 2(m v + m c ) ^xt^max) ( ng ) 
2 fo 2 K max 
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The prefactor f(x) can be taken to be slowly varying and its value at x nmx can be used. 
Integrate the exponentials using the gaussian integral e~ x2 dx = ^pn to obtain the final 
expression for the WKB tunneling probability, 



abs,cm/ rn\ I kg 



AQ\ML \ 2 (m v m r 



1 2 



exp (-2 f*™J dx'K c (x'] £)) 
ZT^/^x-E) 



exp (2 ^7 dx'K v (x>; E ± hu; ko ) / 

£-dx'/K c (x';E±hw ]i0 ) {JlU] 
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